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Accreting X-Ray Binaries display a wide range of behaviours. Some of them 

■ are observed to spin up steadily, others to alternate between spin-up and spin- 
down states, sometimes superimposed on a longer trend of either spin up or spin 
down. Here we interpret this rich phenomenology within a new model of the 
disk-magnetosphere interaction. Our model, based on the simplest version of 

CN ■ a purely material torque, accounts for the fact that, when a neutron star is in 

■ the propeller regime, a fraction of the ejected material does not receive enough 
energy to completely unbind, and hence falls back into the disk. We show that the 
presence of this feedback mass component causes the occurrence of multiple states 

o available to the system, for a given, constant value of the mass accretion rate M* 

q ■ from the companion star. If the angle \ of the magnetic dipole axis with respect 

■ to the perpendicular to the disk is larger than a critical value Xcrit, the system 
eventually settles in a cycle of spin-up/spin-down transitions for a constant value 
of M* and independent of the initial conditions. No external perturbations are 

^ ■ required to induce the torque reversals. The transition from spin up to spin 

down is often accompanied by a large drop in luminosity. The frequency range 
spanned in each cycle and the timescale for torque reversals depend on M*, the 
magnetic field of the star, the magnetic colatitude x, and the degree of elasticity 
regulating the magnetosphere-disk interaction. The critical angle Xcrit ranges 
from ~ 25° — 30° for a completely elastic interaction to ~ 40° — 45° for a totally 
anelastic one. For \ ^ Xcrit, cycles are no longer possible and the long-term 
evolution of the system is a pure spin up. We specifically illustrate our model in 
the cases of the X-ray binaries GX 1+4 and 4U 1626-67. 

Subject headings: accretion, accretion disks — binaries: close — stars: neutron 
— stars: magnetic fields — pulsars: general 
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1. Introduction 

Accreting X-ray binaries, with luminosities up to ~ 10 38 — 10 39 erg/s, constitute the 
brightest X-ray sources in the sky. Since their discovery (Giacconi et al. 1971) more than 
three decades ago, they have provided a unique laboratory to study, among other things, the 
physical processes regulating accretion onto a strongly magnetized neutron star (typically, 
B > 10 11 G). In neutron star (NS) binary systems containing a supergiant or a low-mass star, 
mass transfer takes place through Roche Lobe overflow, and the specific angular momentum 
of matter is sufficiently high to form an accretion disk (the conditions for forming a disk are 
however not necessarily met in Be-star systems, where the NS accretes from capturing the 
star's wind; e.g. Rappaport 1982; Henrichs 1983). In disk-fed systems that host a strongly 
magnetic neutron star, the material from the disk is channelled toward the magnetic poles, 
where it releases its gravitational energy giving rise to the X-ray luminosity we observe. 
Pulsations at the neutron star spin frequency are thus generated from a lighthouse-like 
effect. While the luminosity yields an estimate of the mass accretion rate, pulse timing 
measurements allow one to measure the torque, and hence probe the nature of the accretion 
process mediated by the magnetosphere of the star. Early works (Pringle & Rees 1972; 
Davidson & Ostriker 1973; Lamb, Pethick & Pines 1973) showed that, when accretion occurs 
through a prograde disk, the angular momentum transferred by the accreting material to 
the star (material torque) tends to spin the star up, until the centrifugal barrier inside the 
corotation radius of the magnetosphere (Illianorov & Sunyaev 1975) becomes large enough 
to inhibit further accretion. The star is then expected to settle in a state with an equilibrium 
spin period which depends mainly on the mass accretion rate provided by the companion 
and the neutron star magnetic field (e.g. Frank, King & Raine 2003). 

Observations of disk-accreting X-ray pulsars during the 1970s and 1980s were rather 
sparse, and appeared to be roughly compatible with the near-equilibrium picture (e.g. Nagase 
1989), although there were already hints at times of some unexpected behaviours. These 
included torque reversals for some time while still accreting, or spin up rates much smaller 
than expected for the observed luminosity. In the 1990s, continuous monitoring of several 
disk-fed X-ray pulsars with the Burst and Transient Source Experiment (BATSE) on board 
the Compton Gamma-Ray observatory, shed light on the long-term behaviour of several 
objects (see Bildsten et al. 1997 for a comprehensive review). Particularly striking were the 
findings for GX 1+4 and 4U 1626-67: after about 15 years (for GX 1+4) and 20 years (for 
4U 1626-67) of spin up, both systems showed a torque reversal, which made them switch 
to a spin-down phase. Other systems, like Cen X-3, Vela X-l, Her X-l, often showed an 
alternation of spin up and down sometime overimposed on a longer term of either spin down 
or spin up. In most cases, the magnitude of the torques is comparable during the spin-up 
and the spin-down regimes. These unusual behaviours were a sign that the simple scenario 
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outlined above might be incomplete, and hence they triggered a revival of research, mostly 
in the direction of finding other sources of torque in addition to the one provided by the 
accreting material alone. 

Gosh & Lamb (1979a,b, GL) and Wang (1987, 1995) suggested that, in addition to the 
material torque, there is also an extra source of torque provided by the magnetic field lines 
threading the disk. While in the model of Pringle & Rees (1972) the disk is truncated at 
the point at which the magnetic pressure of the magnetosphere balances the pressure of the 
accreting material, in the GL model there exists a broad transition zone in which the magnetic 
field lines still thread the disk even if the viscous stress in the disk material dominates over 
the magnetic stress. This is made possible through the combination of a number of effects, 
such as the Kelvin-Helmotz instability, turbulent diffusion and reconnection. 

In the models of Arons et al. (1984) and Lovelace et al. (1995) the extra torque is 
provided by the expulsion of a magnetically-driven wind. Transitions between spin up and 
spin down states are possible, but they must be induced by external perturbations, such 
as variations in the viscosity parameter a of the disk or, most plausibly, the accretion rate 
from the companion star. These variations would have to be finely tuned just so that the 
two torque states have comparable magnitude but opposite sign. This seems unlikely in 
general, but even more so in a system like 4U 1626-67, in which the average mass accretion 
rate is likely determined by the loss of orbital angular momentum via gravitational radiation 
(Chakrabarty et al. 1997a). Alternatively, in the case of GX 1+4, Makishima et al. (1988) 
and Dotani et al. (1989) suggested that the spin down could be due to accretion from a 
retrograde disk formed from the stellar wind of the red giant companion. White (1988) 
however showed that this was unlikely to be the case. A retrograde disk around the NS 
spin axis could also be produced by magnetic torques generated in the interaction between 
surface currents on the disk and the component of the NS magnetic field parallel to the disk 
(Lai 1999). 

In this paper we discuss a new scenario for the spin up/spin down transitions observed 
in binary systems accreting from a disk. The torque exchange between the magnetosphere 
and the disk material is supposed to be dominated by the material component as in the early 
models (Pringle & Rees 1972). In this respect, our toy model is very simple and idealized: 
possible torques non parallel to the rotation axis are neglected, as well as magnetic torques 
(e.g. Gosh & Lamb 1979; Lai 2003). What is new in our model is a computation of the 
fate of the ejected material during the propeller phase of the neutron star. Our calculation 
accounts for the following facts: i) not all the "propelled" material receives sufficient energy 
to unbind from the system; ii) if the magnetic moment of the neutron star is inclined with 
respect to its rotation axis, there can be, at the same time, regions of the magnetospheric 
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boundary which are allowed to accrete while others are propelling material away. This is 
a fundamental assumption of our model. While in this paper we provide arguments in 
its support, a final validation will have to wait for detailed numerical simulations. This 
work should therefore be considered as an investigation (the first of its kind to the best 
of our knowledge) of the characteristic timing behaviour of a pulsar whose magnetosphere 
can simultaneously eject and accrete matter in different regions of its boundary. As it will 
be shown in the following, accounting for this possibility leads to fundamentally different 
conclusions for the long-term, equilibrium state of the system. Rather than settling at the 
equilibrium period at which the Keplerian frequency of the disk matches the star rotation 
frequency at the point of interaction (e.g. Frank, King & Raine 1997), the system settles, 
for a wide range of conditions, in alternating cycles of spin-up/spin-down for a constant 
accretion rate from the companion star. A qualitative summary of our model is described 
below, and is formalized mathematically in the following sections. 

A magnetic neutron star surrounded by an accretion disk is able to accrete only under the 
condition that the velocity of the magnetosphere at the point of interaction (magnetospheric 
radius, Rm) is smaller than the local Keplerian velocity of the disk material. If this condition 
is not satisfied, accretion is inhibited (Illiaronov & Sunyaev 1979), and angular momentum 
is transferred from the star to the gas. Whether this propelled gas can be completely 
unbound from the system will depend on the location of the magnetospheric radius within the 
gravitational field of the neutron star. There exists a minimum distance, i? in f, beyond which 
ejection of matter to infinity is possible. If R M < i? in f, the propelled material cannot be 
unbound, and therefore it will fall back on the disk and accrete again. This matter is, in this 
sense, "recycled" . An accreting system with recycled material can, under certain conditions, 
have multiple states available. This is due to the fact that, for the system to be in a steady- 
state condition, the total mass inflow rate at the magnetospheric boundary (which determines 
the position of the magnetospheric boundary itself), M tot = M acc + M rec + M e j e must be such 
that M acc + M e je = M*, where M* is the mass inflow rate provided by the companion star, 
and M acc , M rec and M e j e are, respectively, the rate at which mass is accreted, recycled and 
ejected. Whenever the term M rec is non- negligible, there could be in principle different 
solutions to the above condition corresponding to the same value M* of the mass inflow rate. 
As the system spins up or down on a certain branch of the solution, this solution can be 
lost, and the system is consequently forced to jump to a different state, often characterized 
by opposite torque. This qualitative argument is formalized mathematically in detail in §2, 
while §3 presents specific applications to the cases of the accreting sources GX 1+4 and 
4U 1626-67. Our results are summarized and discussed in 54. 
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2. Model description 

2.1. Magnetosphere-disk interaction in an oblique rotator 

In this section we discuss the main concepts and assumptions upon which our disk- 
magnetosphere model is based. The basic geometry is depicted in Figure 1. The axis of the 
magnetic moment p of the neutron star (NS) is inclined with respect to the rotation axis 
by the magnetic colatitude x- I n cylindrical coordinates (r,<f),z), where the z-axis coincides 
with the rotation axis, the component of the magnetic field in the disk plane is ( Jetzer et al. 
1998) 

5 2 = ^[l + 3(sin X sin0) 2 ], (1) 

under the assumption that the disk is planar and its axis is parallel to the spin axis of the 
NS. When the rotation axis of the NS is inclined (i.e. X 0)> the strength of the magnetic 
field in the plane of the disk depends on the longitude 0. As shown below, this angular 
dependence results in an asymmetric magnetospheric boundary. 

The fate of the matter funnelled from the accretion disk to the rotating, magnetized 
neutron star depends on a number of factors, the most important of which are the relative 
strength of the magnetic pressure and the pressure of the accreting material, and the relative 
velocity of the magnetosphere of the star at the inner radius of the disk with respect to 
the Keplerian velocity at that same radius. Following Lamb & Pethick (1974), the former 
condition can be formalized by equating the magnetic energy with the kinetic energy of the 
infalling matter: 

1 2 B 2 

= 85F ■ < 2 > 

In the free-fall approssimation the density is given by p = pff = M / (4tiv f /r 2 ) , where 
Vff = (2GM/r) 1 / 2 is the free-fall velocity. Using these expressions, together with Eq.(l) and 
(2), the magnetospheric radius for an oblique rotator can be obtained (Jetzer et al. 1998; 
see also Campana et al. 2001): 

R M ((f)) = 3.2 x 10 8 p% 7 M~ 1/7 M~ 2/7 [1 + 3(sin x sin 0) 2 ] 2/7 , (3) 

where p^Q is the magnetic moment in units of 10 30 G cm 3 , Mi is the NS mass in units of 1 M & 
and M17 is the accretion rate in units of 10 17 g s^ 1 . The minimum radius -Rm(O) also cor- 
responds to x — 0, the approximation usually adopted in models of the disk-magnetosphere 
interaction. The maximum radius Rm(^/2) is only a factor of (1 + 3 sin 2 x) 2 ^ 7 < 1-49 larger. 
Note that the elongated shape of the magnetospheric boundary plays a fundamental role in 
our model. 
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An important assumption of our model is that, during the rotation of the magnetosphere 
(whose shape depends upon the instantaneous position of the magnetospheric radius as a 
function of 0) , matter in the Keplerian disk is able to fill the region that separates the disk and 
the magnetospheric flow on a timescale shorter than the spin period of the star. This ensures 
that the inner boundary of the disk remains in constant contact with the magnetosphere. 
We show in the appendix that the Kelvin-Helmholtz instability operates on a sufficiently 
short timescale and wide range of radii that this assumption can be justified. 

Accretion to the star is possible only under the condition that, at the magnetospheric 
radius, the Keplerian velocity of the accreting gas, £Ik(Rm), is larger than the velocity 
Q of the rotating magnetosphere of the star (equal to the velocity of the star), otherwise 
centrifugal forces will inhibit accretion (Illianorov & Sunyaev 1979). The above condition 
is equivalent to saying that the magnetospheric radius must be smaller than the corotation 
radius, -R co = (GM/fip) 1 / 3 , which is the radius at which the Keplerian frequency of the 
orbiting matter is equal to the NS spin frequency Q . In an oblique rotator, the onset of the 
propeller stage will occur when Rm(4>) = R C o at least in one point of the magnetospheric 
boundary. Note that, while a parallel rotator can be either in the propeller or in the accreting 
regime, an oblique rotator can be in both states simultaneously for different longitudes of 
the magnetospheric boundary. Indeed, this special feature of the oblique rotator was used 
by Campana et al. (2001) in building up a model that explained the dramatic luminosity 
variations seen in the BeppoSAX observation of the transient X-ray pulsar 4U 0115+63 1 . 

The interaction between the magnetosphere of the NS and the matter in the disk is 
likely to be at least partially anelastic because of dissipative effects in the mixing process 
between the magnetospheric plasma and the disk matter during the propeller phase. For 
clarity of presentation, here we consider first the two limiting cases of a completely anelastic 
and a completely elastic interaction, and then we will generalize our results to the partially 
anelastic case. 

In the anelastic case, the magnetic field of the NS is able to force matter to corotate at 
the same velocity of the star, and it is endowed at the magnetospheric boundary with specific 
kinetic energy e = 1/2VLqR 2 m and angular momentum / = VLqR 2 m . In order for matter to be 
ejected from the system via the propeller mechanism, the magnetic field must provide it with 
enough energy to reach a velocity in excess of the local escape velocity at Rm- Because in 
the anelastic case the ejection velocity is v e j = QqRm, the requirement above converts to an 



1 This simultaneous presence of different regimes, which is crucial to our model, has not yet been seen 
in numerical simulations. However, to the best of our knowledge, current numerical simulations of the 
propeller regime (e.g. Romanova et al. 2004) are axisymmetric; because of this geometry, they cannot verify 
the simultaneous presence of different regimes of the kind discussed here. 
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" ejection radius" 

i?inf,ane = (2GM/S)*) 1 ' 3 ~ 1.26i? co . (4) 

Only matter which is located beyond this radius during the interaction with the magne- 
tosphere of the NS can be unbound from the system through the propeller mechanism. 
Therefore there exists a region (R co < Rm < Rmf) in which the propeller is active but mat- 
ter cannot be unbound from the sytem by merging with the disk matter (Spruit & Taam, 
1993). We assume that matter in this zone is swung out and circularizes at the radius where 
its angular momentum equals the Keplerian value, i.e. when I = VLqR 2 m — Ik — ^k{Rk)Rk 
(here Rk is the circularization radius). This condition defines the Keplerian circularization 
radius: 

i?K,ane - • (5) 

Matter that is not ejected from the system will fall back into the disk and restart its motion 
toward the NS from the radius defined in Equation (5). 

In the case of the elastic propeller, we assume that material in the disk at Rm moves 
toward the magnetosphere with a tangential relative velocity of — v xe \ = Rm{^o — ^k), 
where fix is the Keplerian angular velocity at Rm- In a completely elastic interaction this 
matter bounces off at the magnetospheric boundary with an opposite velocity of +i> re i that 
in the non-rotating frame sums with t> rot = VLqRm- Thus the ejection velocity is v e j = 
Rm [2fi — Q K (R M )]- In this case the requirement that this velocity be larger than v esc (RM) 
can be written as: 
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where we have used the definition of the corotation radius. This equation can be solved as 
a function of the magnetospheric radius to define the limit beyond which ejection of matter 
to infinity is possible in the purely elastic case: 



R 



inf ,el 



1 + \/2l 2 / 3 

' " R co ~ 1.13^ . (6) 



The matter leaving the magnetospheric boundary is endowed with specific angular momen- 
tum l e i = R M (2Q — Qk)', equating this to Ik gives a new circularization radius for matter 
that is not ejected to infinity. Using the same notation as above we find: 

^k,el - m • (7) 

Let us now consider the most general case of a partially elastic interaction. Following the 
formalism developed by Eksi et al. (2005), we define the "elasticity parameter" (3, which is 
a measure of how efficiently the kinetic energy of the neutron star is converted into kinetic 
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energy of ejected matter through the magnetosphere-disk interaction. Taking into account 
the definitions given above, we now consider the generalized rotational velocity of matter at 
the magnetospheric boundary: 



where 7 = (1 + /3)(1 — Qq/Qk)- The elastic case is obtained in the limit (3=1, and the totally 
anelastic one when /3 — 0. Using Equation (8) we can then generalize also the expression for 
the infinity radius 



In our model we will consider the general case of a partially elastic interaction, and use (3 as 
one of the model parameters. 

Figure 1 (b) shows the various characteristic radii defined above on the disk plane z — 0. 
Depending on the phase (0) and the inclination angle (x), it is possible to have regions of 
the magnetospheric boundary in which accretion is possible {Rm{4>, x) < Rco) together with 
other portions in which the propeller is already active, resulting in ejection of matter to 
larger radii (R co < Rm(<P, x) < Rinf), or to infinity (Rm > Rmf)- In those cases in which the 
inclination angle is sufficiently large, it is possible to have all the three regimes described 
above simultaneously. 

It should be noted that, in our model, we consider ejection of matter from regions of 
the disk that are away from the corotation radius, where the Keplerian velocity of matter 
becomes rapidly supersonic (e.g. Frank, King, Raine 2003). This could in principle lead to 
the formation of supersonic shocks which can heat the plasma and eventually stop the ejection 
mechanism. However in this situation, due to the high relative rotation rate between the 
plasma inside the magnetosphere and that inside the disk, the Kelvin-Helmholtz instability 
can be very efficient. As previously discussed, this instability can lead to a large mixing of 
the two fluids, providing a mechanism to mantain the interaction between the magnetic field 
of the NS and the matter in the disk. Under these circumstances, it has been shown that 
outflowing bubbles of matter are likely to be accelerated magnetically by the NS towards the 
outer region of the disk (Wang & Robertson 1985), in turn supporting the idea that ejection 
far away from the corotation radius can be sustained. 




(8) 




(9) 



and for the circularization radius 



R K = R M (l-l? ■ 



(10) 
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2.2. Conditions for the existence of a limit cycle 

Let M* be the rate of inflowing matter, regulated through the Roche Lobe overflow or 
capture of part of the wind of the companion star. We assume that this matter possesses in 
all cases enough angular momentum that a prograde accretion disk forms. We further assume 
that the mass inflow at the inner disk boundary is azimuthally symmetric (i.e. independent 
of 0). As illustrated in Figure 1 and discussed in §2.1, for a general, oblique, orientation of 
the magnetic field of the NS with respect to the normal to the disk and the spin axis of the 
NS (which we assume are parallel), there will be regions where Rm{4>) < R co , and therefore 
some matter is able to accrete, regions for which Rm{4>) > Rmf that result in matter being 
ejected, and intermediate zones with R co < Rm{4>) < Rmf from which matter gets "recycled". 
The fraction of material in each of these regions is expected to be proportional to the angle 
subtended by the relevant region in the magnetosphere as shown in Figure 1. As in §1, let 
us define M acc , M e j e and M rcc to be respectively the rates of accreting, ejected and recycled 
material at any given time. These various components are illustrated in Figure 2. If dM tot /d(j) 
is the total rate of matter exchanged at the magnetosphere-disk boundary per unit angle, 
these components are given by M comp = 1/27T d<\> {d Mtot/ dcj^j , where the integration 
interval [0i,02] of is such that /?m(0) < R co when comp="acc" , R co < /?m(0) < -^inf 
when comp="rec" and Rm(<P) > Rmf when comp="eje". Figure 3 shows an example of 
these components as a function of the total mass inflow across the entire magnetospheric 
boundary, M tot . At low values of M tot , i? in f > Rm for all values of 0, and therefore all matter 
is ejected (i.e. M e j e = M tot ). On the other hand, at high values of M tot , R co > Rm for any 
0, and therefore all matter is accreted (M acc = M tot ). For values of M tot such that -Rm(0) 
crosses R co at some values of 0, M rcc ^ 0. 

While the total mass inflow rate available to the system is determined by the mass 
transfer rate from the companion, M*, the value of the magnetospheric radius Rm, on 
the other hand, is determined by the total pressure of the accreting matter, i.e. M tot = 
M acc + M e je + M rec . Since, in general, M tot > M*, the magnetospheric radius can be smaller 
than it would be if the "recycled" mass component were not accounted for (as commonly 
assumed in the literature). Therefore, including M rcc in the computation of Rm, allows 
accretion at the same rate to occur for smaller values of M* than it would otherwise. 

In order to demonstrate the existence of a limit cycle, testified by a hysteresis-like loop 
in the M — M tot plane, we start by noting that the rate at which matter is "recycled", M rec , 
does not contribute to the mass budget; therefore a steady-state solution is possible only if 

M* = M acc + M e je • (11) 

Let us therefore examine the behaviour of the curve M tot as a function of the accretion rate 
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M acc + M e je = M*- An example of such a curve for a rotator inclined by an angle of x — 50° 
is shown in Figure 4. All the characteristic parameters of the NS (B, Q , Rns, M ns ) and 
the angle x are kept fixed while the mass supply from the companion is varied. For a given 
value of the external rate of mass supply M*, the corresponding points on the curve yield the 
value (or values) of M tot for which there exists a solution. Again, we stress that the "state" 
of the system, and the characteristics of the solution, are determined by M tot since it is this 
quantity (and not M*) which determines the position of Rm- There can be multiple solutions 
for a given M*, and the one that is realized at a certain time depends on the previous history 
of the system. This situation is reminiscent of a system with hysteresis, and in fact, as 
Figure 4 shows, the shape of the curve M to t(M») resembles a hysteresis curve, where the 
role of the external magnetic field is played by the rate of mass supply by the companion, 
M* (the independent variable in the present context). If at a certain point the system is in, 
say, the state indicated by the point "C" in the figure, and M* increases, the solution (i.e. 
only available state for the system) will be forced to jump to the state indicated by point 
"D" . As M* decreases, the solution will move from "D" to "A" but from that point on, any 
further decrease in M* will cause the solution to jump to point "B". Therefore, like in the 
traditional hysteresis cycle, continuous variations in M* result in discontinuous states for the 
system. 

In the following section, after discussing the computation of the torque, it will be shown 
that the points where the solution jumps from one place to another in the M tot (M#) curve 
often straddle the point of torque reversal. Therefore, transitions between different states 
are often characterized by a torque reversal. 

The case we have illustrated in Figure 4 is only an example of a cyclic behaviour. The 
shape of the curve M tot (M*) changes with the parameters x an d P (while v and B only 
cause a translation in the M tot — M* plane). This can result in several types of cycles with 
a different number of jumps. More examples are shown in Figure 5. 



2.3. Torque and luminosity in the different states of an oblique rotator 

We calculate here the net specific angular momentum trasferred between the disk and 
the NS. In the region of the magnetospheric boundary where accretion is allowed, the net 
specific angular moment transferred from the disk to the NS is given by 

/acc = ^/ (GMi? M ) 1/2 #. (12) 

Z1T JRm<Rco 

In the ejection region, the NS accelerates the material to the ejection velocity, which, as 
discussed in §2.1, is different in the two limiting cases of a completely elastic or anelastic 



- 11 - 



propeller. For the general case of a partially elastic interaction, using Equation (8), the 
angular momentum given by the NS to the ejected matter is: 

ZlX JR M >Rco 



n K ^(l + /5)(fio/fiK-l)#. (13) 

Rm >Rco 

By relating this transfer of angular momentum between the NS and the disk to the variation 
of the NS angular momentum, we have 



~dT ~ 7 



(14) 



where / to t is the sum of the angular momentum computed from Equations (12) and (13), 
and we have assumed that the variation of the NS moment of inertia (/) is negligible. Using 
equations (12) and (13) in (14), we obtain: 

r/L>„ M, 



'tot 



/*Z7T 

/ (GMR M ) l/2 {l - 9(R M - R co )[(l + P)(n /n K - 1)]}# (15) 
Jo 



dt 2tt 

where 9(Rm — R co ) is 1 for R M > R co and for R M < R co . 

Next we compute the different contributions to the luminosity. A schematic represen- 
tation of these contributions is shown in Figure 7. Let us consider first the region of the 
magnetosphere in which there is accretion (Rm{4>) < R co )- The accretion luminosity is given 
by the potential and kinetic energy released by matter falling from the magnetospheric radius 
to the surface of the neutron star; this is 



La.cc 



R M <R C< 



GM (J--±^ + in 2 (i& - R% s ) 



dM.cc , (16) 



where M acc is the fraction of M to t which accretes. Next we consider the contribution to 
the luminosity coming from the "recycled matter". This can be calculated by summing the 
luminosity derived from the release of energy of matter impacting the disk at Rk, and the 
luminosity released from the same matter spiralling in the disk from R K back to Rm- This 
gives 

jR co <R M <Ri I ,f \ z Zfi «/ 
where M rec is the rate corresponding to the "recycling" part of the magnetospheric boundary 

(Rco < Rm(4>) < -Rinf)- 
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Another contribution to the total luminosity is provided by the release of energy in 
the boundary layer which separates the magnetosphere from the Keplerian disk. This term 
applies to matter at any longitude if we consider a completely anelastic propeller ((3 = 0), 
because in this case the magnetosphere forces matter to corotate with it during both the 
accretion and the propeller regime. On the other hand, in the limit of a completely elastic 
propeller (f3 = 1), this term is present only for those angles for which Rm < R co and 
matter is thus slowed down in the boundary layer before it can begin falling toward the NS. 
If we consider a generic value for the elasticity parameter, the luminosity of the boundary 
layer can be written as: 

/ [R 2 M (n 2 K - nl)]d<P for R M < R co 

Lbl = < £ J ° ^ (18) 



m r 

-^(1-/3) J [R M (^o ~ «k)]# for Rm > Rco ■ 



47T 

Finally, we have to account for the luminosity produced by the matter inflowing from the 
companion as it spirals in towards the magnetospheric radius in the Keplerian disk. This 
contribution, which is obviously present in all different regimes, is given by 

It is important to emphasize that in our model both the torque and the luminosity 
depend on the total mass inflow rate M tot at the magnetospheric boundary, and this can 
take different values for the same mass accretion rate M*. The three panels of Figure 5 show 
the behaviour of the torque and luminosity as a function of M tot for three combinations of 
NS parameters. These are chosen to represent different types of limit cycles (also shown 
in the figure for each case - note the axes here are swapped with respect to Figure 4 for 
consistency with the other panels). In Figure 5(a), a transition between the points A and B 
is accompanied by a reversal from spin up to spin down, while the jump from point C to D 
will cause a transition from spin down to spin up. The luminosity is at its lowest at point 
B and at its highest at point D, but the overall variation during the cycle is well within an 
order of magnitude. A more complicated cycle is depicted in Figure 5(b); here a transition 
from point A to B causes a spin-up to spin-down reversal, while the opposite happens during 
the jump from point E to F. This cycle comprises also another jump, from point C to D, 
with both points on the spin-down branch. The luminosity varies by more than three orders 
of magnitude during the cycle, being at its lowest during most of the spin down phase. The 
third example of limit cycle, the one shown in Figure 5(c), has only two allowed jumps, 
both of them straddling the point of torque reversal, as in case a), but the luminosity is 
substantially larger when the system is on the spin- up branch (A - D), than when it is on 
the spin-down branch (B - C). 
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Whether there exists a limit cycle depends crucially on the angle x- this has to be 
large enough to ensure that some regions of the magnetosphere are in the accretion regime 
while, at the same time, others are in the propeller phase. There exists a critical value of 
the magnetic colatitude, Xcrit, below which the steady state solution breaks into two disjoint 
curves and it is no longer possible to find a cyclic behaviour through a sequence of steady- 
state solutions. This is illustrated in Figure 6 for the cases (3 — and (5 — 1. If the accretion 
rate M* from the companion is above a certain value (which depends on u, B, x and (3), only 
one solution is available to the system, and it corresponds to the spin up branch (see Figure 
5). On the other hand, if M* < M crit (for the example under consideration, M CT it = 2.4 x 10 16 
g/s for (3 = and M crit = 3.5 x 10 16 g/s for (3 = 0, but it varies with v and B), then multiple 
solutions are available for any value of x displayed, and the one that is realized at any given 
time depends on the history of the system. However, a cyclic jump of the solutions between 
the spin- up and the spin-down branches can only be realized for angles above Xcrit- 

The critical angle ranges from about 25° — 30° for (3 — 1 to about 40° — 45° for (3 = 0, 
and, for a given (3, it is independent of v and B. Therefore, for the curves shown in the figure, 
a limit cycle can only be achieved in the cases with \ — 30° and x = 40° for (3 = 1, and in 
the cases with x — 50° and x — 70° for (3 = 0. In the other cases displayed, the M*(M tot ) 
curve is discontinuous. The shape of the curve is such that, if the system is spinning up, a 
solution on the spin- up branch can be found for any value of M*, and therefore the system 
will continue spinning up. On the other hand, if the system is originally on the spin-down 
branch (which is possible only for < M crit ), then any decrease in will keep the system 
on the spin-down branch, while an increase in above M CT it will cause a jump on the 
spin-up branch, and from that point on the system will be spinning up independent of the 
value of M*. Note that, depending on the angle x-, there can be spin-up solutions even at 
very low mass inflow rates M*. This result is a novelty of our model, deriving from the fact 
that the recycled mass component M rec can keep the magnetospheric radius in "pressure" 
even if M* is very small. 

Among all the components that make up the total luminosity, the accretion term is the 
only that is certainly pulsed, since the accretion material is funnelled by the magnetic field of 
the star onto the NS magnetic poles, where its energy is released. The accretion luminosity 
therefore varies with the phase of the star, resulting in a pulsating flux. Also the boundary 
layer luminosity might be pulsed at the NS spin. Therefore, the maximum pulsed fraction 
in our model is constrained to be between / pul = L acc /L tot and / pu i s = (L acc + L Bh )/L tot . 

Note that the sum of the various contributions in Equations (16), (17), (18), (19) can 
result in a complex, non-monotonic dependence of L tot as function of the accretion rate 
from the companion, M*. In the classical model of accretion onto magnetized neutron stars, 
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the transition between the standard accretion regime onto the NS surface to the regime 
of accretion onto the magnetospheric boundary in the propeller regime is marked by the 
change between the oc M* and the oc M^ 7 scaling of the luminosity (Stella et al. 1994; 
Campana & Stella 2000). In the propeller phase, the underlying assumption of these works 
is that the main contribution to the luminosity derives from the disk luminosity (Eq. 19). 
In the present model, this might not be the case if there is a non-negligible contribution 
to the luminosity from recycled matter. Moreover, at low accretion rates, we find that the 
contribution to the luminosity from the boundary layer generally dominates over that from 
the disk for an anelastic propeller (see Figure 8). For sufficiently low values of M* (so that 
^k(Rm)/^o <^ 1), the boundary layer luminosity scales as oc M* 3/7 , while L BL oc Ml' 7 
at high values of M* (for which Vt^(R M )/^% 3> 1). When the corotation radius is of the 
order of the magnetospheric radius, however, these dependences are changed. Since the 
Keplerian frequency at the magnetospheric radius is an increasing function of M*, in the 
propeller regime the term |J1q — Q^(Rm)\ decreases with the increase of M*, while in the 
accretion regime the same term increases with increasing M*. As a result, when Rm is of 
the order of R co , the luminosity of the boundary layer has a flatter dependence on M* for 
R M > R co and a steeper dependence for R M < R co . This can be seen in Figure 8. Both 
the case of a completely anelastic propeller ((3 = 0), and a totally elastic one ((3 = 1) are 
considered, showing respectively the maximum and the minimum boundary layer luminosity 
that the system can have. In the former case we find that, for sufficiently low accretion 
rates (so that the whole magnetospheric boundary is in the propeller regime), the boundary 
layer luminosity is substantially larger than the disk luminosity. The relative contribution 
LBiV-^disk dearly increases as the degree of anelasticity increases, since L BL oc (1 — j3). As 
the mass accretion rate M* increases, so that at least some regions of the magnetospheric 
boundary are in the accretion regime, the disk luminosity begins to dominate over that of 
the boundary layer (this is now independent of (3). However, L BL has a stronger dependence 
on M*, and, for sufficiently large M* that Rm <C R co , L bl becomes ~ -^disk- 

The two panels in Figure 8 show the cases of a slow pulsar [y = 9 mHz) and a fast one 
[y = 100 mHz). The discussion above regarding the relative contribution of Lbl and Ldisk 
to the total luminosity budget holds in both cases. Furthermore, once accretion sets in, the 
accretion luminosity dominates over both Ldisk and Lrl- The slower the pulsar, the larger 
is this term compared to the others. Therefore, in the accretion regime and for Rm -C R co , 
L tot oc M* as in the classical models. However, for Rm ~ -Rco, the presence of the "recycled" 
term of luminosity in our model causes a non-monotonic dependence of the total luminosity 
on M*, with multiple solutions allowed. The actual solution that is realized at any given time 
will depend on the history, i.e. whether the system is on the spin-up or spin-down branch of 
the limit cycle (see Figure 4). This is an important difference of our model with respect to 
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the classical solution where, once the system is in the accreting phase, the luminosity scales 
monotonically with M*. On the other hand, Figure 8 shows that there are regions for which 
a small variation in M* can cause a large jump in luminosity. 

Similar to the luminosity, the behaviour of the torque in the surroundings of the region 
with Rm ~ -R co is complex and non-monotonic. Small variations in M can cause the system 
to jump between states with opposite sign of the torque. Within this region, because of the 
complex dependence of both L and Cl on M*, our model does not make any specific prediction 
regarding correlations between torque and luminosity. In most situations, these are expected 
to be uncorrelated, and different types of limit cycles (see Figure 5) will generally lead to 
different behaviours in the various spin-up and spin-down phases. 

2.4. Cyclic spin-up/spin-down evolution at a constant M* 

In systems in which mass transfer takes place through Roche Lobe overflow, the rate 
at which material is fed to the disk is expected to be roughly constant, or characterized by 
relatively low-amplitude, long-term variations. We are not concerned in this section with 
the accretion disk instabilities that likely give rise to the very large amplitude variations of 
the mass inflow rate in binary X-ray transient systems. Rather, in the following we describe 
how recurrent episodes of spin up and spin down can be achieved in our model in response 
to a strictly constant accretion rate M* from the companion star. 

Figure 9 shows the behaviour of the curve M = M e je + M acc (top panel) and the cor- 
responding frequency derivative, z>, (bottom panel) as a function of M tot and for different 
values of the period (corresponding to different times). The parameters B, x and (3 are the 
same in all cases. They yield a limit cycle of the type described in Figure 5(b). While the 
specific points of torque reversal will vary depending on the type of cycle (as shown in the 
various examples of Figure 5), the underlying structure determining the transitions is the 
same in all cases and therefore we analyze in detail only one of the possible scenarios. 

In order to illustrate how the spin up/spin down states are achieved at a constant M*, 
let's start, say, with the system at a frequency v so that the corresponding M(M tot ) curve is 
the one labeled "2" in Figure 9, and let's assume that the system is in a spin- up state. The 
intersection between the curves M and M* on the spin-up branch of the cuspid determines 
the value of Mtot, M tot sol , corresponding to the allowed spin- up state for that value of the 
frequency. This value of M totiSO i in turn determines the value of the frequency derivative at 
that point in time (point "O" in both panels of the figure). The frequency at time t + dt is 
simply determined as vit + di) = u{t) + dz/(M totso i)/<it. As the pulsar spins up, the curve 
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"2" moves towards curve "1" until the point at which the spin-up branch of the solution 
rises higher than the system M* (point A). From that point on, the only possible state for 
the system that satisfies the condition M = M* is the one corresponding to point C in the 
figure, on the spin-down branch (negative torque). Once again, the new (current) value of 
A^tot.soi determines the actual value of v (corresponding point C in bottom panel) which is 
used for the next time step to determine the new v. While on the spin down branch of the 
solution, the curve M now moves from the curve "1" towards the curve "2" and then "3" . 
Spin down continues until this branch of the solution does not intersect any longer the M* 
line (point B), at which point the only allowed state for the system to be is on the spin-up 
branch, and the system reverses from spin down to spin up. This is the beginning of a new 
cycle. 

In this model, the points of spin reversals are determined by the maximum and minimum 
of the M(v) curve. The shape of this curve depends on the anelastic parameter (3 and on the 
angle x- For a given (5 and Xi a change in the strength of the magnetic field simply results 
in a shift of the curve without a change in shape: a higher B field would move the curve 
to higher values of M tot , therefore resulting in stronger spin-up and spin-down torques, and 
hence in a shorter timescale for torque reversals. For the model to work as described, it is 
clear that the points where the solution jumps must straddle the point of torque reversal. We 
find this to be the case for a wide range of combinations of x an d j3. However, for each value 
of (3, there is a narrow range of angles x f° r which the torque inversion point falls outside 
the allowed region for the transitions. For these particular and rare cases, the system would 
tend towards the point Cl = and remain there, for a strictly constant M*. However, small 
fluctuations in M* can still cause the system to jump from one solution to another. For 
the rest of this discussion we will focus on the greatest majority of cases for which torque 
reversals naturally occur at M* = const, unless we explicitly state otherwise. 

If the magnetic colatitude angle x is larger than Xcrit, the system is bound to end up in 
a cyclic sequence of spin-up/spin-down transitions. In fact, as it can be seen from Figure 9, 
if the system starts with a much larger frequency than the maximum frequency in the cycle, 
i/ max , it will spin down since only one solution (on the spin down branch) is allowed as long as 
v > f max . Similarly, if the system starts with a frequency much smaller than the mimimum 
frequency in the cycle, iVin, it will spin up as only one solution (on the spin up branch) is 
allowed as long as v < z/ min . Therefore, our model predicts that the system, independent of 
the initial conditions, eventually settles in a region where there are cyclic transitions between 
spin-up and spin-down states. This limit cycle is not induced by external perturbations, but 
is the natural equilibium state forwards which the system tends. 
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3. Application of our model to persistent X-ray pulsars 

In the following, we will apply our model to two objects for which long term- monitoring 
showed a marked transition between a spin-up and spin-down phase. We will then discuss 
the way our model can be generalized to other cases where short-term episodes of spin- 
up/spin-down are superimposed onto longer term spin-up or spin-down trends. The most 
comphrensive monitoring of the spin behaviour of accreting X-ray pulsars in binaries is given 
in Bildsten et al. (1997), and here we briefly summarize the observations for the two cases 
that we model. 



3.1. GX 1+4 

GX 1+4, discovered in 1970 through an X-ray balloon experiment (Lewin, Ricker & 
McClintock 1971) is an accreting X-ray pulsar binary hosting an M red giant (Davidsen et 
al. 1977); the orbital period is likely to be of a few years (Chakrabarty & Roche 1997). 
Early observations through the 1970s showed that this source was spinning up at a very 
high pace with a spin-up timescale \vji>\ ~ 40 yr. The frequency changed from ~ 7.5 mHz 
to ~ 9 mHz during the first 15 years of observations. In the early 1980s, however, the flux 
dropped abruptly and the source could not be detected by Ginga. Given the sensitivity of 
the instrument, the flux must have decresead by more than two orders of magnitudes for a 
few years. Once its flux raised, the source could be monitored again, and it was found to 
spin down on a timescale comparable to the previous spin-up timescale (Makishima et al. 
1988). 

A solution that closely reproduces the observed source behaviour was found by running 
the time-dependent code described in §2.4 for a range of parameters B, x, P- The correspond- 
ing value of M* is determined so that the point of torque reversal of the system between spin 
up and spin down matches the observed value. The larger the magnetic field, the larger M* 
and hence the torque, and therefore the more rapid the timescales of the torque transitions 
are; the parameters x an d (3, by determining the shape of the M curve, especially influence 
the total frequency range z/ max — u min that the system spans in a cycle. 

For the case of GX 1+4, we found that a good choice of parameters is the combination 
B = 6 x 10 13 G, x — 45°, and (3 = 0.3. These yield a cycle of the type displayed in panel (b) 
of Figure 5 and in Figure 9. In particular, the parameters B, x and (3 used in Figure 9 are the 
same as those used for GX 1+4. The accretion rate provided by the donor companion must 
be M* ~ 2 x 10 16 g/s in order to produce a turnover in frequency around 9 mHz. With this 
choice of parameters, Figure 10 shows the behaviour of the system that our model predicts. 
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Cycles of spin- up/spin-down alternate in response to torque reversals. The luminosity of the 
source is comparable during the spin-up and spin-down phases, except for a few years around 
the time of spin reversal from spin up to spin down, when it drops abruptly. This is due to 
the fact that, after the system has "jumped" to point "C" in Figure 9 (at the beginning of the 
spin down phase), there are no regions in the magnetosphere-disk boundary where accretion 
onto the NS can take place, the NS is in the propeller regime, and therefore M acc = (see 
Figure 3). During that time, the only contribution to the luminosity comes from the disk and 
the boundary layer, which are however much smaller than the accretion luminosity (since 
this is a slow pulsar). A prediction of our model is that, while large drops in luminosity 
can be expected when the system reverses from spin up to spin down, they should not 
occur in correspondence of the spin-down/spin- up transition, because when this transition 
occurs (refer to the jump from point "B" to point "A" in Figure 9), most regions at the 
magnetospheric boundary are allowed to accrete. While these overall features are generally 
robust predictions of our model, the detailed variation of M (and hence the luminosity) 
with torque shown in our examples should not be taken too rigorously. These variations 
depend on the shape of the M(M tot ) curve, and this is in turn determined by the shape 
of the magnetospheric boundary as a function of time. As discussed in §3.3, a number of 
effects neglected here can influence this shape, and hence affect the detailed behaviour of the 
solution. In particular, note that observations of GX 1+4 show that luminosity and torque 
strength are correlated during part of the spin down phase (Chakrabarty et al. 1997b). This 
feature is not reproduced by the current version of our model. 



3.2. 4U 1626-67 

4U 1626-67, discovered by SAS-3 in 1977 (Rappaport et al. 1977) is an ultracompact 
binary with an extremely low-mass companion (Levine et al. 1988; Chakrabarty et al. 1997a) 
and a 42 minute orbital period (Middleditch et al. 1981). During the first ~ 20 years of 
observations, the source was found to spin up with a timescale vjv of about 5000 yr. The 
frequency increased from 130.2 mHz to about 130.5 mHz, at which point the source started 
to spin down. Unlike the case of GX 1+4, there was no evidence for a large change in the 
bolometric luminosity of the source during the transition. 

The very long timescale for spin reversal of this source (due to a smaller torque compared 
to the case of GX 1+4) requires a smaller magnetic field. We found that our model yields 
a reasonable match to the observations with the choice of parameters B = 2.5 x 10 12 G, 
X = 68°, (3 = 0. The corresponding solution found with our model is displayed in Figure 11. 
The upper panel shows only one spin- up/spin-down torque reversal, since the complete spin- 
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up/spin-down cycle, of the order of several thousand years, lasts much longer than the 
observed time. Although the luminosity somewhat drops around the time of spin reversal, 
it does so to a lesser extent and for a much shorter time than for the case of GX 1+4. The 
reason for these differences lies in the variation of the shape of the function M(M tot ) for 
different choices of the parameters \ and (3. The parameters that best match the solution 
for 4U 1626-67 yield a cycle of the type in panel (a) of Figure 5. The transition from spin 
up to spin down (point A to B in the figure) is accompanied by a less dramatic variation 
in luminosity than it is for the cycles of the type shown in panels (b) and (c). Note how, 
for this source, since the observation window is much smaller than the timescale for torque 
reversal, other torque inversions are not expected in the near future, unless induced by 
external perturbations. 



3.3. Generalizations and limitations of our model 

The two examples given above, for two sources spinning up and down at very different 
rates, show that our model can reproduce different types of cyclic behaviours. In the two 
cases discussed, we assumed that the mass accretion rate from the companion, M*, does 
not vary with time. Under this assumption, our model predicts that the points of torque 
reversals will always occur at the same value of the frequency. On the other hand, if the 
donor accretion rate varies with time, this will no longer be the case. If M* increases 
with time, then the points of torque reversals will occur at larger frequencies as time goes 
on. Viceversa if M* decreases with time, then the points of torque reversals will occur at 
smaller and smaller frequencies with time. A combination of discrete states in an oblique 
rotator (producing cyclic torque reversals), with longer-term variation in the external M* 
can produce a long-term spin evolution with superimposed shorter cyclic episodes of spin up 
and spin down. 

Also note that, depending on the system parameters (namely the inclination angle x 
and the elasticity parameter (3, which determine the shape of the M(M tot ) curve, and hence 
the points of torque reversals), the transition from a state of spin up to a state of spin down 
can result in a period of time during which accretion is completely inhibited (i.e. M acc = 0) 
and the luminosity is orders of magnitude lower (unless the pulsar has a very fast spin in the 
ms range and the luminosity of the disk and the boundary layer are conspicuous even when 
M acc = 0). The system can then behave as a "transient" even when the accretion rate from 
the companion is constant. 

In the present (simplest) version of our model, the frequency range (kWx - ^min) spanned 
in a spin-up/down cycle cannot however be made arbitrarily small. In order for the torque 
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reversals to occur at constant M* and without any other external perturbation, the curves 
M(z/ max ) and M(z/ min ) (curves 1 and 3 respectively in Figure 9) must be such that the two 
points of torque reversals (A and B in Figure 9) satisfy the conditions MA(f max ) = M* and 
MB^min) — respectively. Arbitrarily small cycles require arbitrarily small loops in the 
M(Mtot) curve, so that the inversion points can be extremely close. This cannot be achieved 
with the current version of our model, in which the shape of the M(M tot ) curve (and hence 
the "size" of the loop around the points of torque reversals) depends only on the inclination 
angle x an d the elasticity parameter (5. However, there are a number of effects that we have 
neglected here, and which could be potentially important for small-scale torque reversals. In 
particular, if the disk plane is not orthogonal to the NS rotation axis, a precession of the disk 
around the spin axis can be induced (Lai 1999), producing a time- dependent modulation of 
the various regimes on a time scale on the order of the spin period of the star. We reserve 
to future work a more comprehensive exploration of the physical effects that influence the 
magnitude and frequency of the torque reversals. 



4. Summary and Discussion 

A magnetic rotating neutron star surrounded by an accretion disk is an intuitive example 
of an accreting system in which the conditions can be realized such that a fraction of the 
matter is accreted, another fraction is ejected and completely unbound from the system, and 
another part is propelled out but does not possess enough energy to unbind, and therefore 
falls back onto the disk, getting "recycled". We have shown that, for a given mass rate 
supply from the companion, accretion with the mass feedback term included leads to multiple 
available states for the system, characterized by different (and discrete) values of the total 
mass inflow at the magnetospheric boundary. The luminosity in each of these states is 
generally different, as it depends on the relative amounts of the various components of the 
total mass inflow rate. The available states often straddle the point of torque reversal, and 
therefore correspond to states with opposite sign of the torque. 

The character of the solutions is essentially determined by the inclination angle x of the 
NS axis with respect to the disk. At angles x ^ Xcrit, the limit cycle breaks down. In this 
case, for an external mass supply larger than a critical value (which depends on the system 
parameters), the system can only be on the spin up branch. For accretion rates smaller than 
this critical value, both the spin up and the spin down branches of the solution are possible, 
and the one that is realized will depend on the history of the system. After a sufficiently 
long time, however, if the system is spinning down, the available solutions will be drifting 
and the source will jump out of the spin-down branch, and continue evolving on the spin- up 
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branch. For x ^ Xcrit, cyclic transitions between states of opposite torque can be realized 
even at a constant value of the accretion rate from the companion. This is a particularly nice 
feature of our model: periodic variations between spin up and spin down states take place 
without requiring the presence of any external, periodic, and fine-tuned perturbation. Most 
importantly, we have shown that periodic, cyclic episodes of spin up/spin down behaviour 
must be realized in a number of situations. While in the classical theory of accreting X-ray 
binaries (where the effect of mass feedback is not accounted for) the system is expected 
to eventually settle at the equilibrium frequency which matches the Keplerian frequency at 
the magnetospheric boundary, in our model, where recycling is accounted for, the system 
will eventually settle around a limit-cycle behaviour in which different spin derivative and 
luminosity states alternate, recurrently. The points of spin reversal and the timescales of 
the torque reversals depend on a combination of factors, namely the accretion rate from the 
companion, the magnetic field of the NS, the inclination angle of the NS axis, and the degree 
of anelasticity at the disk-magnetospheric boundary. 

In the case of the two X-ray binaries GX 1+4 and 4U 1626-67, we have determined a set 
of parameters B,x, (3 which is able to reproduce the main features of their timing behaviours, 
such as the timescales and frequency span of the transitions, as well as the large luminosity 
drop observed around the transition spin up-down in the case of GX 1+4 but not of 4U 1626- 
67. The correlation between torque strength and luminosity in the spin down phase observed 
in GX 1+4 (Chakrabarty et al. 1997b) is however not reproduced by the present scenario. 
On the other hand, we still need to emphasize that ours is a very simplified model and 
therefore the detailed behaviour of our solution should not be considered too rigorously: 
while our model appropriately accounts for the material torque at the disk-magnetospheric 
boundary when a fraction of mass is recycled, it neglects other possible sources of torque, 
such as magnetic stresses (e.g. GL) or magnetically driven outflows in an extended boundary 
layer (Arons et al. 1984; Lovelace et al. 1995). The presence of other torque terms could 
modify the character of the solutions if non-material torques dominate over the material one. 
A general treatment that includes all possible sources of torques is beyond the scope of this 
paper, especially since the relative strength of the various terms would be hard to estimate 
from first principles. 

Finally, while the details of the solutions that we have discussed specifically apply to 
the case of a rotating neutron star accreting from a disk fueled by a companion star, the 
general feature of a multeplicity of states available for a given mass inflow rate of matter 
can probably be generalized to other accreting systems in which "recycling" occurs. An 
example is that of an accretion disk around a rotating black hole. Numerical simulations 
(e.g. Krolik et al. 2005) show that, while a fraction of the accreting mass is ejected through 
a jet, another fraction, of slower velocity and at larger angles from the jet axis, falls back 
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into the disk, getting recycled. It would be interesting to include this mass feedback process 
into numerical simulations of accretion disks around black holes, and investigate whether 
the discountinuos states and cyclic behaviour might ensue in those cases as well. 
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APPENDIX 

Here we justify our assumption that, during the rotation of the magnetosphere, matter 
in the disk is able to fill the region that separates the disk and the magnetospheric flow on 
a timescale shorter than (or comparable to) the spin period of the star. 

Let r v = R 2 /u = R/vr be the viscous timescale in the disk, where R is the radial 
distance from the star, v the kinematic viscosity coefficient and vr the radial velocity in the 
disk. In the reference frame of the disk (in which r v is measured), the stellar rotation time is 
T rot = 27r/jQo — S^k] - Using the thin disk approximation, the disk height H can be written as 
H = f'R, where / << 1 is a numerical factor that can be assumed approssimatively costant 
for small variations of the radial distance from the NS (typically / ~ 1/10). Furthermore, 
using the a prescription for the viscosity (Shakura & Sunyaev 1973), we can write vr = 
ctVg/vK = af 2 VK where v s is the sound speed in the disk and we have assumed v s /v K ~ H/R. 

Let us consider first the propeller regime (Rm > R co )- For R = Rm, we obtain that 
t v < r rot only if 

R co <R M <{1 + 2ixaf 2 ) 2 '*R co (1) 

which corresponds to a very narrow region around the corotation radius. Beyond this region, 
the viscous timescale becomes too long to permit a replenishment of the inner regions of the 
disk as the star rotates. Here another mechanism is needed to justify our assumption. Indeed, 
in the propeller regime, the surface of separation between the magnetospheric and disk flow 
is Kelvin-Helmholtz (KHI) unstable due to the large shear velocity (Wang & Robertson 
1985; Spruit & Taam 1993). In the frame corotating with the NS this velocity is v re \ = 
Rm[Q* — 0, k (Rm)]- Because of the KHI, matter in the disk is mixed with the NS magnetic 
field lines, thus mantaining a strong interaction between the disk and the magnetosphere. 

The characteristic timescale for the development of the KHI (in the direction of the 
shear motion) can be estimated as (e.g. Stella & Rosner 1984) r KH ~ 47r(/c|v re i|)~ 1 , where 
k = 2n/\ is the wave vector of the perturbation which inizializes the instability. The 
condition that the KHI develops within a time shorter than the local timescale r rot /2 is 
hence satisfied for wave vectors k > 2/R. Furthermore, in order for the interaction between 
the disk and the magnetospheric flow to be maintained throughout the rotation of the star, 
the KHI must be able to mix disk matter and magnetic field lines at least on a distance 
d ~ [R M (0) - Rm{k/2)} < 0.5R M (0) (see Eq.(3) and Fig.l). The simulations of Wang 
& Robertson (1985) show that perturbations of lengthscale A become rapidly unstable and 
evolve into elongated vortices of magnitude comparable to A. This means that a perturbation 
of length A is able to produce mixing between matter and field lines on a distance scale of the 
same order. Wang & Robertson also argue that the dominant mode of the instability will 
likely be the one just sufficient to offset the effect of viscous damping through the turbulent 
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motions in the shear layer. In our case this condition traslates into \/2nv t ~ (f2* — Qk)^ 1 
where v t is the turbulent velocity. If we choose v t ~ v s and use / ~ H/R, we can roughly 
estimate 

that is of the order required to cover the radial extension d discussed above (here we have 
used the fact that in our model the propeller regime typically occurs for Rm within the 
range R co < R M < (1.6 — 1.7)i? co ). Since the wavelength in Eq(2) satisfies the condition 
k > 2/R M , the dominant mode of instability develops in a shorter time than the local 
dynamical timescale, and therefore the KHI is able to maintain a close interaction between 
the disk and the magnetosphere on this timescale. 

Let us now consider the accretion regime (Rm < R co ). Using the same derivation as 
above, the analogous of Eq.(l) is 

(1 - 2iraf 2 f 3 R co < R M < R co (3) 

which is again a narrow region in the vicinity of the corotation radius. In the accretion 
regime however, considering the argument used in the propeller case (where now t> re i = 
Rm(^k — ^o))j we obtain the same conclusion about the efficiency of the KHI, and the 
analogous of equation 2 is now 

which clearly satisfies the requirement A > 0.5Rm for any value of Rm in the region of 
interest. Furthermore, after the KHI has brought matter just inside the magnetospheric 
radius, the enhanced contribution of the gravitational with respect to the centrifugal force, 
forces matter to fall toward the NS also under the effect of the Rayleigh- Taylor instability 
(Arons & Lea 1980; Wang & Robertson 1984). This enhances the transport of matter toward 
the NS and therefore strengthens the reliability of our assumption. 
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Fig. 1. — Schematic illustration of the NS-disk system for an oblique rotator. Figure (a) 
shows the relative positions of the magnetic dipole moment axis, the phase angle 0, and 
the inclination angle \. The NS is assumed to rotate around the z-axis. Figure (b) is a 
two dimensional rapresentation of the position of the magnetosphere (continuous line) with 
respect to the corotation radius (long dashed line) and the infinity radius (short dashed line), 
for arbitrarily fixed values of the NS parameters. 
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Fig. 2. — The fate of the matter provided by the companion at a rate M* depends on the 
relative position of the magnetospheric radius with respect to the corotation radius and the 
ejection radius. Matter can be accreted, ejected or recycled into the disk. 
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Fig. 3. — The various contributions to the total accretion rate M tot = M acc + M rec + M e - ]e at 
the magnetospheric-disk boundary. The system parameters are B = 6 x 10 13 G, v — 9 mHz, 
X = 45°, (3 = 0.3. 
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Fig. 4. — Schematic representation of the hysteresis-type limit cycle. The arrows indicate 
the points where the system "jumps" between different states as a result of variations in 
the external mass supply rate M*. The system parameters in this example are v — 9 mHz, 
B = 6 x 10 13 G, x = 80°, p = 0. 
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Fig. 5. — Various types of hysteresis limit cycles. The system parameters are v — 9 mHz, 
B = 6 x 10 13 G and (3 = in all cases, while x = 80° i n panel (a), \ = 50° in panel (b) and 
X = 47° in panel (c). In the top panels of each case, the arrows indicate the points where the 
system "jumps" between different states as a result of variations in the external mass supply 
rate M*. The middle panels show that, under most circumstances, a jump is accompanied 
by a torque reversal and, in some cases, by an abrupt change in luminosity (displayed in the 
bottom panels). 
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Fig. 6. — The break of the cyclic solution in the M* — M tot plane is shown for a system 
with v — 9 mHz s, B = 6 x 10 13 G, /3 = (left panel) and (3 = 1 (right panel). When the 
inclination angle is small, it is no longer possible to find a steady state, cyclic solution. The 
value of x around which the solution breaks depends on the anelasticity parameter (5 but is 
independent of the values of v and B. 
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Fig. 7. — The various contributions to the total luminosity budget for an accreting neutron 
star. 
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Fig. 8. — The various contributions to the total luminosity budget L tot = L acc + L TCC + Ldisk + 
Lbl for an accreting neutron star as function of the mass accretion rate from the companion, 
M*. The system parameters are v = 9 mHz, B = Q x 10 13 G, x — 80° in the top panels and 
v — 9 mHz, B = 10 9 G, x — 80°; m the bottom ones. The two limiting cases of a completely 
elastic interaction (left panels) and of a completely anelastic interaction (right panels) are 
shown. 
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Fig. 9. — Variation of the function M(M tot ) = M acc + M e j C (top panel) and of the correspond- 
ing spin rate variation (bottom panel) at three different times during a spin up/spin-down 
cycle. The system parameters are .8 = 6x10 13 G,x = 45° and (3 = 0.3. The three curves cor- 
respond to frequencies v = 9.1 mHz (curve 1), v = 8.3 mHz (curve 2) and v = 7.7 mHz (curve 
3). Allowed states for the system are only those satisfying the condition M* = M(M tot ); 
in this case, log(M*) = 16.3. When multiple solutions are allowed, the state in which the 
system will be found depends on its previous history. 
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Fig. 10. — An oblique NS rotator with magnetic field B = 6 x 10 13 G, inclination angle 
X = 45 cleg, and elasticity parameter f3 = 0.3 is able to reproduce the main spin-up/spin 
down characteristics of GX 1+4. The luminosity is comparable during the spin- up and spin- 
down phases, except for a few years at the beginning of the spin-down phase, when it drops 
abruptly. 



-37- 




Fig. 11. — An oblique NS rotator with magnetic field B = 2.5 x 10 12 G, inclination angle 
X = 68°, and elasticity parameter j3 = is able to reproduce the main spin- up/spin down 
characteristics of 4U 1626. 



